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Resistive jet simulations extending radially self-similar 
magnetohydrodynamic models 
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ABSTRACT 

Numerical simulations with self-similar initial and boundary conditions provide a link 
between theoretical and numerical investigations of jet dynamics. We perform axisym- 
metric resistive magnetohydrodynamic (MHD) simulations for a generalised solution 
of the Blandford & Payne type, and compare them with the corresponding analyti- 
cal and numerical ideal-MHD solutions. We disentangle the effects of the numerical 
and physical diffusivity. The latter could occur in outflows above an accretion disk, 
being transferred from the underlying disk into the disk corona by MHD turbulence 
(anomalous turbulent diffusivity), or as a result of ambipolar diffusion in partially 
ionized flows. We conclude that while the classical magnetic Reynolds number i? m 
measures the importance of resistive effects in the induction equation, a new intro- 
duced number, Rp = (j3/2)R m with /3 the plasma beta, measures the importance 
of the resistive effects in the energy equation. Thus, in magnetised jets with (3 < 2, 
when Rp < f resistive effects are non-negligible and affect mostly the energy equa- 
tion. The presented simulations indeed show that for a range of magnetic diffusivities 
corresponding to Rp > 1 the flow remains close to the ideal-MHD self-similar solution. 
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1 INTRODUCTION 

Collimated outflows of plasma observed to emerge from the 
vicinity of a wide spectrum of cosmic objects are still a chal- 
lenge for observational and theoretical astrophysics. These 
outflows play a key role in the transport of angular momen- 
tum and energy of the accreted gas facilitating thus, for ex- 
ample, star formation. Nevertheless, when new observations 
put more and more severe constraints on models, these seem 
to be still too rudimentary to provide sophisticated answers. 

The starting point of the modeling of jets are the ideal 
MHD equations, which can be solved analytically by assum- 
ing axisymmetry, time-independence and the self-similarity 
ansatz. Analytical models of ideal MHD disk winds (Bland- 
ford & Payne 1982, re-visited in Vlahakis et al. 2000; here- 
after V00) , provide not only the first insight into the physics 
of such outflows but equally important they can be used as 
a test bed of more sophisticated simulations of the resis- 



* E-Mail: mikiOtiara. sinica. edu.tw (MC); jgracia@cp.dias.ie 
(JG); vlahakis@phys.uoa.gr (NV); tsingan@phys.uoa.gr (KT) 



tive MHD system via various numerical codes. In Vlahakis 
& Tsinganos (1998) general classes of self-consistent ideal- 
MHD solutions have been constructed. Two sets of exact 
MHD outflow models have been found: meridionally and ra- 
dially self-similar ones. Previously known studies were recog- 
nised to belong in this more general classification of all avail- 
able analytical models. In particular, the V00 study reme- 
died the physically unacceptable feature of the Blandford & 
Payne (1982) terminal wind solution which was not causally 
disconnected from the disk (see also Ferreira & Casse (2004) 
when a resistive disk is included). 

Among the basic problems which still remained to be 
solved however were the common deficiency that all radially 
self-similar models had, namely, a cut-off of the solution at 
small cylindrical radii and also at some finite height above 
the disk where they were unphysical. The reason for such 
behaviour is a strong Lorentz force close to the system's 
axis. The invalid analytical solution very close to the axis 
has been corrected numerically. Also, a search in the numer- 
ical simulations for solutions at larger distances from the 
disk has been performed, in Gracia et al. (2006; hereafter 
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Figure 1. Initial conditions for our numerical simulations. The 
solid lines represent logarithmically spaced isocontours of den- 
sity. It is also shown in colour grading. The dashed lines show 
poloidal ficldlines or flowlines. The dotted lines indicate, from top 
to bottom, the position of the fast-magnetosonic {small dots), the 
Alfvcn, and the slow-magnetosonic surface (large dots), respec- 
tively. 



GVT06), where ideal-MHD numerical simulations with the 
V00 solution as initial condition have been performed using 
the NIRVANA code (version 2.0, Ziegler 1998). These results 
have been verified also by using the PLUTO code (Mignone 
et al. 2007) in the ideal MHD simulations by Matsakos et 
al. (2008). 

The next step in the exploitation of the available ana- 
lytical solutions has been the investigation of the dynamical 
connection of the outflow to the underlying disk (see e.g., 
Konigl, 1989; Wardle & Konigl 1993; Li 1995; Ferreira 1997; 
Casse & Keppens 2004; Zanni et al. 2007), providing some 
understanding of the formation of jets from the accretion 
disk. At the same time, since the disk is naturally resistive, 
the need emerged to go beyond the ideal MHD regime. With 
the magnetic field included, the effects of the magnetic re- 
sistivity, both numerical and physical, had to be addressed, 
discriminated and analysed. 

In fully ionized disks an anomalous turbulent diffusivity 
has to be invoked to allow for accretion of matter crossing a 



large-scale magnetic field. This anomalous turbulent diffu- 
sivity may be present in the outflow as well, at least at dis- 
tances close to the disk. In partially ionized disks the physi- 
cal conductivity is properly described by a tensor, consisting 
of three distinct parts corresponding to the ambipolar diffu- 
sion, the Hall effect, and the Ohmic dissipation (e.g. Wardle 
& Ng, 1999; Salmeron et al., 2007). The dominant mecha- 
nism in the outflow above the disk is most likely the ambipo- 
lar diffusion (e.g., Sano & Stone 2002; Kunz & Balbus 2004; 
Wardle 2007), and can be appropriately described only by 
a multifluid MHD. Nevertheless, in both cases (turbulent or 
ambipolar diffusion), a scalar conductivity can capture the 
basic characteristics of the breakdown of ideal MHD related 
to the magnetic field diffusion and the resistive heating. 

The effects of the resistive heating in the formation and 
acceleration of jets from resistive disks or tori, have been 
studied in Kuwabara et al. (2000, 2005), concluding that 
Joule heating is not playing an essential role in jet forma- 
tion. However, in these studies only low values (lower than 
the critical value that we define below) of resistivity were ex- 
amined. Resistive effects have been also studied in Fendt & 
Cemeljic (2002; hereafter FC02). In this case, however, the 
energy equation was not solved and a polytropic equation 
of state has been assumed instead, such that the effects of 
finite resistivity have not been directly incorporated in the 
energetics of the problem. Safier (1993a) working in the am- 
bipolar diffusion regime of outflows associated with young 
stars, found that the diffusive term in the energy equation 
cannot be neglected. The heating of the gas could have sig- 
nificant observational consequences (see e.g. Safier 1993b; 
Martin 1996b; Cabrit et al., 1999; Garcia et al., 2001a,b; 
O'Brien et al, 2003; Shang et al, 2004). Safier's (1993a) 
work also shows why the cooling term in the energy equa- 
tion and the ionization balance equation need to be consid- 
ered in a full investigation. He was able to uncover a strong 
feedback mechanism between the gas temperature and the 
ionization fraction (which scales inversely with the ambipo- 
lar diffusion heating rate). The result of the heating depends 
also on the geometry of the flow which strongly affects the 
adiabatic cooling. In a spherical outflow this cooling effec- 
tively counters the Joule dissipation heating (Ruden et al., 
1999), while in a disk-driven jet with a small streamline di- 
vergence the adiabatic cooling is relatively unimportant (at 
least initially). Interestingly, Joule dissipation can play a 
role (although generally not a dominant one) also in mag- 
netically guided accretion problems (e.g., Martin 1996a). 

Numerical resistivity is implicitly present in any nu- 
merical simulation, and its various effects need to be iden- 
tified and studied in detail. This is one of the main aims 
of this paper. The other one is to investigate the effects of 
small physical resistivity, examine the stability of the resis- 
tive jet solutions and define when the resistivity is "small" 
and when it becomes "large". The novel approach followed 
here is that the analytical solution is a stationary reference 
solution, something which was lacking in the previous inves- 
tigations of resistive MHD flows, e.g. FC02. 

The structure of the paper is as follows. The analytical 
expressions and their modification in the setup for numerical 
simulations are first presented in Sec. 2. Then the resistive- 
MHD solutions (Sec. 4) are systematically compared to the 
ideal-MHD ones of Sec. 3. In Sec. 5 we introduce an exten- 
sion of the magnetic Reynolds number which quantifies the 
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transition from an ideal-like behaviour of the low diffusivity 
solutions for values of the diffusivity r\ below a critical value 
r/crit, to a transient and erratic behaviour of the solutions 
for r\ > ri cr it. The main possible consequences for astrophys- 
ical outflows are briefly discussed. A summary of the main 
results is given in the last Sec. 6. 



2 PROBLEM SETUP 
2.1 Governing equations 

The resistive-MHD equations solved by the NIRVANA code 
are, in SI units: 
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where V is the flow velocity, B is the magnetic field, (p, P) 
are the gas density and pressure, and $ = — QM/r is the 
gravitational potential of the central mass M. The internal 
energy (per unit mass) is related to the pressure and density 
by 

e = r - > (6) 

7-1 p 

where 7 is the effective polytropic index. The magnetic diffu- 
sivity r\ is assumed constant and is related to the resistivity 
pc = /io?7, where /10 is the permeability of vacuum. 

2.2 Initial and boundary conditions 

For the initial and boundary conditions of the simulations 
the self-similar solution of V00 is used. The assumptions of 
steady-state, axisymmetry and radial self-similarity result 
in the following expressions for the physical quantities in 
spherical (r, 8, <f>) and cylindrical (Z — rcos8, R = rsin#, 
<h) coordinates: 
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Here we decomposed vector quantities in poloidal (index p) 
and toroidal (index <f>) components. We set the solution pa- 
rameters to (x , A 2 , fi , k , 7) = (0.75 , 136.9 , 2.99 , 2 , 1.05), as 
in the V00 solution. 

The diffusivity 77 (which is assumed constant through- 
out the domain) is normalised as 77 = rj V0R0 = 
fi\/GMRo/n, with 77 dimensionless. The self-similar solution 
breaks down near the rotation axis. This becomes evident 
from the fact that all physical quantities are proportional 
to a power of the function 1/a, which is divergent on the 
axis (R = 0), see Eqs. (|7I12|I . In addition, the analytical so- 
lution of V00 is not provided for 8 smaller than 0.025 rad, 
measured from the axis. 

To perform numerical simulations in a computational 
box with the symmetry axis included, we need to mod- 
ify/extrapolate the analytical solution. Near the axis, we 
extrapolated the missing analytical solutions for the tab- 
ulated functions G, M and ip, as described in GVT06. A 
similar result can be achieved with less involved extrapola- 
tion, as shown in Matsakos et al. (2008). Modification of the 
functions G and M means also that the pressure/energy is 
modified near the axis. 

For the magnetic field in the vicinity of the symmetry 
axis there is an additional problem. With the extrapolated 
functions G and ip the magnetic field as given by Eq. @ is 
not divergence-free. This leads to the need for suitable mod- 
ification of the initial magnetic field. A simple modification 
is to compute the Bz component from the Z component of 
the self-similar expression 
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and subsequently the radial component Br by solving the 
V ■ B — with boundary condition Br(R = 0) = 0. 

The poloidal velocity field should, in result, also be 
modified. Initially V p || B p , as demanded for steady ideal- 
MHD flow. Therefore, we compute the new poloidal velocity 
field, maintaining the magnitude of velocity, but correcting 
the direction so that V v \\ B v holds. 

The modified initial magnetic field is presented in Fig.[T] 
together with the modified density isocontours and positions 
of the critical surfaces. 

For boundary conditions we use symmetry conditions 
on the rotation axis, and outflow conditions on the outer 
R and Z boundaries. On the lower boundary Z = Z m i n 
(we take Z m i n = 6R0 in all simulations) we fix the values 
for six physical quantities, namely density, three velocity 
components, the azimuthal and one poloidal magnetic field 
component (the other one is given from V • B = 0). The 
pressure/internal energy boundary condition is kept fixed 
only in the super-sonic region (where the poloidal velocity is 
larger than the speed of sound), otherwise it is extrapolated 
from the flow onto the boundary. 

More specifically, in our simulations the Z-component 
of the magnetic field in the first ghost cell is set from the an- 
alytical solution, and the radial component is obtained from 
the divergence-free condition. Therefore, the magnetic flux 
along the boundary is fixed at all times. In the second ghost 
cell the Z-component of the magnetic field is extrapolated 
linearly, to allow for change in the field line shape, when the 
radial component is again obtained from V • B = 0. 
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Figure 2. Ideal-MHD simulations before, during and after the relaxation, for the resolution of R X Z = (128 X 256) grid cells = 
([0, 50] X [6, 100])-Ro. Lines denote thirty logarithmically spaced isocontours of density. The times of simulations are few thousands, few 
ten thousands and larger than few hundred thousands of the Courant time steps, Left to the Right panel, respectively. 



In our computations we used various resolutions and 
sizes of the computational domain. Here we present the 
results for the resolution R x Z — (256 x 512) grid cells 
= ([0,50] x [6, 100])7?o, in the uniform grid. Results com- 
ply with the solutions for one fourth, one half and double 
of this resolution, which we also computed and presented, 
when needed for direct comparison. 



to right, the kinetic, enthalpy, Poynting, and gravity terms 
can be recognised. 

The degree of alignment of the lines on the poloidal 
plane where the above quantities are constant together with 
the poloidal magnetic field lines can be used as a test on how 
close to a steady-state is the final result of a simulation. 



2.3 Ideal-MHD integrals 

It can be shown, that steady, axisymmetric, ideal-MHD 
polytropic flows conserve five physical quantities along the 
poloidal magnetic field lines (Tsinganos 1982). These so 
called integrals are the mass-to- magnetic-flux ratio ^a, the 
field angular velocity fi, the total angular momentum-to- 
mass flux ratio L, the entropy Q, and the total energy-to- 
mass flux ratio E (henceforth we call the latter integral en- 
ergy for brevity). These integrals are given as 
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The various contributions of the energy E correspond to the 
various terms on the right hand-side of Eq. (|19[1 . From left 



3 IDEAL-MHD SIMULATIONS AND 
NUMERICAL RESISTIVITY 

Ideal-MHD (rj — 0) numerical simulations with the same 
setup and using the same code have been performed in 
GVT06. The density isocontour plots for ideal MHD sim- 
ulations in different times are presented in Fig. for illus- 
tration of the relaxation process. 

In the early stage of simulation, the initial conditions, 
especially the setup near the symmetry axis, affect the out- 
flow time-evolution. Up to few thousands of Courant time 
steps, the solutions might look moderately different, before 
the relaxation towards the stationary state. Finally, after few 
ten thousand Courant time-steps (or in higher resolutions 
after few hundred thousands) the simulations give similar 
result as the initial state. This solution is stationary (not 
quasi-stationary) as we did not note virtually any change 
in a time ten times longer than the time needed to reach 
the stationary state. This time is equal to five millions of 
Courant times steps, or 2500.Ro/Vb, when expressed in nor- 
malised units. 

The initial setup near the axis, which may introduce big 
differences in the initial state (compared to the self-similar 
solution), does not affect the reached final state so much as 
the corresponding part of the boundary Z = Z m i n near the 
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Figure 3. Illustration of the effect of numerical resistivity, i.e. grid resolution. Left: The slow-magnetosonic, Alfvenic, and fast- 
magnetosonic critical surfaces (from bottom to top). Different line types represent the final states of simulations with resolution 128 X 256 
(dotted/red), 256 X 512 (dot- dashed/blue), and 512 X 1024 (dashed/green). The solid/black lines show the initial-state critical surfaces 
in the high resolution reference simulation. Middle: The shapes of two different magnetic flux surfaces, i.e. poloidal magnetic field lines, 
for the same resolutions and line types as described above. The solid/black lines represent again the initial state of the high resolution 
reference simulation. Right: Same as in the Middle panel, but for the energy (E) integral lines. 
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Figure 4. Illustration of the effect of numerical resistivity, i.e. grid resolution, on the alignment of the MHD integrals with the magnetic 
flux surfaces. Left: As an example for the MHD integrals, the entropy Q, normalised to its value at large distances, is plotted along the 
inner total energy integral line shown in Fig. [3] for the initial state (solid/black) of the high resolution reference solution, and the final 
states at resolutions 128 X 256 (dotted/red), 256 X 512 (dot- dashed/blue), and 512 X 1024 (dashed/green). Right: Same as in the Left panel 
but here Q is plotted along a magnetic field line (instead of the energy integral line). 



disk surface. These ideal MHD results have been extensively 
discussed in GVT06. 

For our purpose here, to use them as a reference solution 
for the resistive runs, the problem of numerical resistivity 
should be addressed. 



In ideal-MHD numerical simulations, a current sheet, 
which forms at any discontinuity of the magnetic field, 
should not diffuse away. Also, magnetic field lines should 
not reconnect through such a sheet. However, since in a 
finite-difference scheme computations can not resolve fea- 
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Figure 5. Illustration of the effect of numerical resistivity, i.e. grid resolution, on the contributions to the energy integral E. Left: 
Split-down of the energy contributions along the two integral lines shown in Fig. [3] The upper set of curves corresponds to the inner 
integral line, the lower set of curves to the outer integral line. Colours code the resolution as in Fig. [3] The different line types represent 
energy E (solid), kinetic energy (dotted), enthalpy (dot-dashed) and Poynting (dashed), respectively. The gravitational energy is not 
shown (it is orders of magnitude smaller). Right Same as Left panel, but plotted along the field line instead of the integral lines. 



tures smaller than grid cells, numerical reconnection occurs, 
cf. Hawley & Stone (1995). 

To check the effects of the numerical resistivity, we com- 
pared our ideal-MHD simulations performed in various res- 
olutions. These were Rx Z = (128 x 256) , (256 x 512) and 
(512 x 1024) grid cells, in identical setups. In Fig.|3]the effect 
of grid resolution (i.e. the effect of the numerical resistivity) 
on the position of the critical MHD surfaces, on the shape of 
the field lines, and on integral lines, is shown. For increasing 
numerical resistivity (i.e. for lower resolution), the critical 
surfaces move downstream. However, towards large cylin- 
drical radii R and small height Z, the effect of the boundary 
conditions becomes important and the critical surfaces bend 
towards the disk, as noted in GVT06. Also, for increasing 
numerical resistivity the field lines tend to straighten out 
and move to smaller cylindrical radius. However these dif- 
ferences are insignificant. 

In Fig.|3]the entropy integral Q = p/p~" is shown, along 
the same field lines as in Fig. [3] normalised to its value at 
large distance. 

In the ideal MHD case, all integral lines should coincide. 
Fig- 51 illustrates that the initial state does not show perfect 
alignment of the integrals. This is expected, since the initial 
state is a modified exact solution of the ideal-MHD equa- 
tions. For the final state the integral lines at any resolution 
are better aligned than in the initial state. As expected, the 
alignment of the integrals is better for lower numerical re- 
sistivity, i.e. higher resolution. In contrast to the analytical 
solution, the nume rical solution features a shock (see a re- 
lated discussion in lMatsakos et al.l 120081 ). As expected, the 
entropy jumps across the shock, as illustrated in Fig. [3] but 
is otherwise constant along the energy integral line. 

Further, for ideal-MHD steady flows the integrals fall 
exactly on magnetic flux surfaces, i.e. field lines. However, 



as shown in the right panel of Fig. [4] the presence of numer- 
ical resistivity changes this situation. While the integrals 
stay well-aligned among themselves, the alignment with the 
magnetic field lines is poor for low resolution (128 x 256) and 
becomes almost perfect for high resolution (512 x 1024) runs. 
In fact, while the integrals re-align while the simulation pro- 
gresses, the field lines diffuse away from the corresponding 
integral lines over time. However, for all numerical resolu- 
tions in our simulations this process eventually comes to halt 
and a stationary state is reached. 

Not only are the integrals aligned very well, but also the 
individual contributions to the energy E along the integral 
lines are very similar across different numerical diffusivities 
as illustrated in Fig. [5] Apart from resolution effects at the 
shock along the inner integral lines, all models follow similar 
trends and converge to the same profiles. 

Again, plotting the same quantities along field lines, as 
shown in the right panel of Fig. [S] shows large spread across 
different numerical resistivities. Higher values of resistivity 
show larger changes along the field line, as well as do field 
lines further in. 

We conclude that numerical resistivity does effect the 
solution, but in a smooth manner. For reasonable resolu- 
tions (grid cells small enough compared to the characteristic 
length of the problem), it does not challenge the solution in 
our setup. 



4 RESISTIVE-MHD SIMULATIONS 

To investigate the resistive-MHD behaviour of such out- 
flows, we set the resistive-MHD numerical simulation for 
NIRVANA with the same initial and boundary conditions 
as for the ideal-MHD simulations presented in the previous 
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Figure 6. Illustration of the effect of physical resistivity. Left: The slow-magnetosonic, Alfvenic, and fast-magnetosonic critical surfaces 
(from bottom to top). Different line types represent the final states of simulations at resolution R X Z = (256 X 512) grid cells, with the 
physical magnetic resistivity i) = {dotted/red), i) = 0.03 {dot-dashed/blue) and i) = 0.15 {dashed/ green) . The solid/black lines show the 
initial state critical surfaces. Middle: The shapes of two different poloidal magnetic field lines, for the same resolutions and line types 
as described above. The solid/black lines represent again the initial state. Right: Same as in the Middle panel, but for the energy (E) 
integral lines. 




Figure 7. Illustration of the effect of physical resistivity on the alignment of MHD integrals and magnetic flux surfaces. Left: As an 
example for the MHD integrals the entropy Q, normalised to its value at large distances, is plotted along the inner energy integral line 
shown in Fig. [6] for the initial state (solid/black), and final states for 17 = (dotted/red), rj = 0.03 (dot- dashed/blue) and r) = 0.15 
(dashed/green) for runs with resolution Rx Z = (256 X 512) grid cells. Right: Same as in the Left panel but plotted along the innermost 
field line (instead of the energy integral line). 



section. The magnetic diffusivity is set to constant through- 
out the computational box. It would be possible to model 
the diffusivity as proportional to a product of a character- 
istic length (e.g. the cylindrical distance) with a character- 
istic velocity of the problem (e.g. the Alfven velocity or the 



sound speed). Alternatively, we could model an ambipolar 
diffusivity as a non-constant scalar 77 = V^t™, where Va 
is the Alfven speed and r„i is the neutral-ion momentum 
exchange time. In the special case that the magnetic field 
and current density are orthogonal (as is the case close to 
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Figure 8. Illustration of the effect of physical resistivity on the contributions to the energy integral E. Left: Split-down of the energy 
contributions along the two integral lines shown in Fig. [5] The upper set of curves corresponds to the inner integral line, the lower set 
of curves to the outer integral line. Colours code the physical magnetic diffusivity as in Fig. \E\ The different line types represent energy 
E (solid), kinetic energy (dotted), enthalpy (dot-dashed) and Poynting flux energy (dashed). The gravitational energy is not shown (it is 
orders of magnitude smaller). Right: Same as Left panel, but plotted along the field line instead of the integral lines. 



the disk where the poloidal field dominates) this is an exact 
expression (e.g., Balbus & Terquem, 2001). However, in this 
study we opted for the simplest case of a constant rj. It is 
expected that, at least for relatively low values of rj, the ex- 
act prescription of diffusivity should not significantly affect 
the flow. As noted in FC02, where diffusivities rj oc p 1 ^ 3 were 
examined!]] differences introduced by relating the resistivity 
to density are small. Variable resistivity cases in which the 
exact prescription may significantly affect the results in the 
high resistivity limit, will be examined in another connec- 
tion. 

The nature of magnetic diffusivity which we introduced 
in our simulations depends on the specific circumstances 
valid for the investigated flow. We treat it as an "effec- 
tive magnetic diffusivity", without discussing its physical 
origin, an issue that would be beyond the scope of this pa- 
per. The most obvious case could be magnetic turbulence, 
which would extend from the disk to the disk corona im- 
mediately above the disk. Since we do not treat the disk, 
but take it as a boundary condition here, we can not treat 
the resistivity self-consistently. Therefore we take it as a free 
parameter. 

We performed a study of magnetic resistivity in our 
simulations. At first, preparatory work has been done, with 
level of numerical magnetic diffusivity tracked by decreasing 
the parameter i) until there was no effect on the solutions 
(i.e. until they became identical to the ideal MHD solutions, 
obtained by the code). In our setup here it showed to be of 
the order of fj ~ 0.001. 



1 These cases correspond to r] oc C s oc p' 7_ oc p 1 / 3 for 
7 = 5/3 and P/ p 1 being a global constant (and not constant 
along the flow as in the self-similar model of V00) . 



Then we started increasing the magnetic diffusivity. The 
solution remained similar in character to the ideal-MHD one 
until some threshold critical magnetic diffusivity fj c = 0.15 
has been reached. For fj > fj c the solution changed abruptly, 
not resembling the initial condition anymore. Similar be- 
haviour has been reported in FC02 who also refer to a crit- 
ical magnetic diffusivity, but in their study a comparison 
to some analytical solution was not possible^ More impor- 
tantly, FC02 ignored the resistive term in the energy Eq. @, 
and thus they could not observe a modification in the flow 
caused by the energy dissipation as we do here. Their crit- 
ical resistivity has to do with the diffusion of the magnetic 
field; as a result we cannot directly compare the two works. 

Analysis and direct comparison of the data for the re- 
sistive runs with the ideal-MHD analytical solutions can be 
performed only when the solutions do not depart largely 
from the stationary ideal-MHD ones. In such case the in- 
tegrals along the similar lines can be compared. For the 
large resistivity, as the solutions differ significantly, a sep- 
arate study of validity and stationarity of the new solution 
is required. 

Therefore, here we concentrate on the solutions not de- 
parting significantly from the ideal-MHD solutions. 

The positions of critical surfaces and field line shapes for 
the different physical resistivities (but for single resolution) 

2 I n FC02 the setup was m otivated by reasons of comparison 
with lOuved fc Pudrit j [119971) ideal MHD numerical simulations, 
which reach well defined quasi-stationary state, needed for suc- 
cessful comparisons, but are not necessarily stationary even in the 
ideal MHD regime. Quasi-stationarity of the solutions in FC02 has 
been defined rather by a "rule of thumb", when here the solution 
reaches well defined stationary state, which is possible to relate 
and compare to the analytical solution. 
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are shown in Fig. [6] As was the case with the numerical 
resistivity, the critical surfaces move downward the flow with 
increasing resistivity. This effect is more prominent here, 
meaning that f) = 0.03 gives an upper limit for the numerical 
resistivity of the R x Z = (256 x 512) resolution. 

Also, the field lines tend to straighten out. However, 
field lines close to the axis seem to be little, if at all, affected 
by the resistivity, leading to the conclusion that the flow is 
not modified along them. This is different than was the case 
for the numerical resistivity. 

The MHD integrals (see Fig. are not well conserved 
along the flux surfaces, as was observed also for the numer- 
ical resistivity. However, the misalignment is not enlarged. 

The evolution of individual contributions of the energy 
(see Fig. [8} shows a clear trend along the flux surfaces. For 
increasing resistivity, the Poynting energy and enthalpy in- 
crease, while the kinetic energy decreases. The differences 
are not big, though, especially in the outer field lines. Again, 
the field lines close to the axis are special in the sense that 
all the curves fall atop of each other, indicating that the 
energy there is independent of the physical resistivity. 

The curves for fj = and fj — 0.03 are almost identical, 
what amounts to the conclusion of our preparatory work 
mentioned above, about the order of numerical magnetic 
diffusivity. It is confirmed now that i) — 0.03 is a low physical 
magnetic diffusivity case. 

As in the numerical diffusivity cases, not only are the 
integrals aligned very well, but also the individual contribu- 
tions to the energy E along the integral lines are very similar 
across different diffusivities, as illustrated in Fig. [8] 

In general, the time evolution of the NIRVANA solu- 
tions for a low magnetic diffusivity parameter fj 0.15 does 
not differ much from the ideal MHD evolution. It only takes 
more computational timesteps to reach the stationary state, 
as the diffusive timestep now adds to the total timestep. 

After a few hundred thousands (or few millions, for 
larger resolutions) Courant timesteps of the relaxation pro- 
cess, the outflow reaches the stationary state, which is sim- 
ilar to the initial state. The relaxation is not dramatic, as 
the flow is without strong shocks, and connection of the ini- 
tial condition to the boundary conditions is smooth. Such 
evolution is expected for given initial conditions. 

Here we also presented solutions for the magnetic dif- 
fusivity close to the highest one which does not change the 
solutions dramatically, fj — 0.15. For larger magnetic dif- 
fusivity than this threshold, the solutions seem to depart 
significantly from the initial condition. Therefore, it is not 
possible to describe these solutions simply comparing the in- 
tegrals along the same lines, as we did in the analysis above, 
because the geometry of the solutions is completely different. 

Preliminary results for the high resistivity regime in our 
setup are shown in Fig. [9] The flow with fj — 1.5 is exam- 
ined, i.e. with diffusivity 50 times larger than the typical 
"low magnetic diffusivity" value. The solution departs sig- 
nificantly from the ideal MHD case, and seems to show some 
periodicity in time-evolution. This result of the radical devi- 
ation of the solution from the stable ideal MHD jet solution 
may have some interesting astrophysical implications. For 
example, one possibility is that if somehow the resistivity 
drastically increases in the disk (e.g., due to some instabil- 
ity), a well formed and behaved jet ceases to exist, giving its 
place to a more erratic outflow. However, this investigation 




50 
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Figure 9. Illustration of the effect of super-critical diffusivity 
on the flow structure. This simulation with fj = 1.5 never reaches 
a steady state, but remains highly time-variable. The solid lines 
represent logarithmically spaced isocontours of density. It is also 
shown in colour grading. The dashed lines show poloidal field- 
lines. The dotted lines indicate, from top to bottom, the position 
of the fast-magnetosonic {small dots), the Alfven, and the slow- 
magnetosonic surface (large dots), respectively. 



is beyond the scope of the present paper, where we check the 
behaviour and stability of the ideal MHD solutions with the 
inclusion of the resistivity for the particular problem setup. 
Instead, it will be examined in a following paper. 



5 A CRITICAL VALUE OF THE DIFFUSIVITY 
IN MHD JETS 

The resistive effects in a magnetohydrodynamic flow 
have been traditionally quantified by using the magnetic 
Reynolds number _R m , the dimensionless ratio of the advec- 
tion and diffusion parts in the induction Eq. )3). In a case 
of a jet the advection velocity is the flow speed V, while a 
characteristic scale measuring the distance at which the var- 
ious quantities vary substantially is the cylindrical radius R. 
Thus, the Reynolds number is, 

VR 

Rm = — , (20) 
V 
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Figure 10. Value of 0.5/3(V R/VoRo) for the analytical solution 
of V00. This quantity gives the critical value r) c that corresponds 
to Rp = 1. 



and magnetic diffusion is unimportant for R m 2> 1. 

However, the resistivity affects also the energy trans- 
port, Eq. Q, through the Joule heating term. The ratio 
of non-resistive terms in this equation over the Joule heat- 
ing term gives another important dimensionless number, 
Rp = (PV/ R) / (r/B 2 / fj,o R 2 ) , or, in terms of the plasma beta 
P = 2ji P/B 2 , 



Rg 



P VR _ P 
2 ~ ~ 2 



Rn 



(21) 



Energy dissipation is important for Rp < 1. Clearly, for 
magnetised jets with p < 2 it is Rp < 7? m . If i? m < 1 
then Rp < 1 as well, meaning that resistivity effects are 
important in both, the induction and the energy equation. 
However, there is a possibility to have i? m 3> 1 and Rp < 1. 
These inequalities define a regime where resistivity affects 
the energy, but not the induction equation. 

Taking flows with progressively larger diffusivity, as the 
numerical experiments that we have performed in this study, 
we will first reach a point where Rp « 1, while still the 
Reynolds number is Rm 3> 1. From this point on, the flows 
will significantly differ from the ideal-MHD initial condi- 
tions, something which is connected to the critical diffusiv- 



(22) 



ity that we observe. By writing r\ = rjRoVo we find that the 
condition Rp = 1 gives a critical value 

Vc ~ 2 V R ' 

This quantity for the analytical solution of V00 is shown in 
Fig. 1101 Inside the region Z/Zo > 6 that we simulate, the 
minimum value of 0.5f3(V R/VoRo) is 0.13 (near the point 
R/Ro ~ 17, Z/Rq = 6). This means that only for r) > 
0.13 resistive effects start to play a role - their influence is 
more important (at least initially) near the point R/Ro « 
17, Z/Rq = 6. This value is very close to the numerically 
evaluated critical value r) c — 0.15. 

Fig. llll shows the behaviour of i? m during the simulation 
with the low magnetic diffusivity f) — 0.03, and Fig. ll2l shows 
Rp. In the captions of these figures listed are also the values 
of the R m and Rp at the bottom of the flow. 

The corresponding minimum values for the simulation 
with super-critical rj shown in Fig. [9] are R m = 5 and Rp = 
0.1 at the bottom of the flow. These values, however, are not 
constant in time since this solution is not stationary, with 
the dense "wing" sweeping the computational box quasi- 
periodically. 

Astrophysical jets, which are presumably launched from 
the accretion disk around a central object, are present in 
various scales and around objects of various masses. The 
magnetic resistivity of the disk - an essential part of the ac- 
cretion mechanism in some models - would be transported to 
the disk corona immediately above the disk. Therefore, the 
effects investigated in this paper can be related to astrophys- 
ical objects. We concentrate on the case of jets associated 
with young stellar objects, but since our numerical simula- 
tions are scalable to objects of any mass, similar scaling to 
the jets around e.g. a black hole, is possible (as long as the 
flow remains non-relativistic) . 

Let us now scale our simulated solutions to the case 
of jets associated with young stellar objects. This is readily 
done by inspecting the equations in the 321 As we have nor- 
malised rj by VoRo, scaling to our physical system is straight- 
forward. 

The velocity Vo is related to the mass of the central 
object by the first of Eqs. (|13|l . This leaves us with the ex- 
pression for rj 

Vgmro 



r) = fjVoRo 



(23) 



Therefore, for any object of mass M we define the unit ra- 
dial distance Ro, and we can estimate the physical rj in our 
simulations. 

For young stellar objects, with mass M « Mq and 
characteristic distance of ~ 0.1AU, i.e. Ro = 20Rq = 1.4 x 
10 m, in our setup, rj = 6.8r) x 10 14 m 2 s _1 . This gives for 
the physical diffusivity, rj = 2 x 10 13 m 2 s _1 for our small r) 
value of 0.03, and rj = 10 14 m 2 s _1 for the last subcritical 
value of r) — 0.15. 



6 SUMMARY 

In this paper we presented resistive numerical simulations of 
outflows with radially self-similar initial conditions. The an- 
alytical solution of V00 has been modified, following GVT06, 
in a way that it provided a consistent setup for simulations 
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Figure 11. The magnetic Reynolds number i? m = VR/n for f\ = 0.03 in the computational box when R X Z = 256 X 512 grid cells, for 
the physical domain R X Z = ([0, 50] X [6, 100])i?o- I n the Left panel the isocontours of R m are shown. The values of R m are from 1500 
to 30000 in 20 contour lines, increasing linearly from bottom (near the disk surface) to top. In the Right panel we show R m for the same 
case, as a function of height above the disk, at a constant cylindrical distance corresponding to the middle of the computational box. 
The minimum R m along this line is 335. 



which aimed to extend the failing analytical solution in the 
close vicinity of the axis, and large distances from the disk. 
These ideal MHD solutions have been confirmed also by 
Matsakos et al. (2008), by using the PLUTO code, which 
is using different numerical methods than the NIRVANA 
code, used in our simulations. 

From the outset, it is not obvious at all that the resis- 
tive MHD solutions for such a problem should stay close to 
the ideal MHD solutions. However, we find that the MHD 
solution changes smoothly, with a continuous trend for the 
physical variables, as the resistivity increases. Resistive solu- 
tions also reach a well defined stationary state. This in itself 
is already an interesting and valuable result of the present 
study. The topology of the moderately resistive MHD solu- 
tions turns out to be similar to the ideal MHD case. 

This is the case until some critical magnetic diffusivity is 
reached, when the solutions become increasingly nonconser- 
vative in energies and fluxes. The critical transition can be 
measured through a new dimensionless quantity Rp, which 
measures the influence of resistive effects in the energy equa- 
tion. 

Energy and flux considerations in the present paper are 
of unprecedented exactness, when it goes for stationarity 
of the compared runs, which can help to reach some conclu- 
sions on resistive MHD behaviour of jets in the astrophysical 
context. Especially it could be of value for the treatment of 
the disk corona nearby the disk, where the resistivity of the 
disk is probably transported up to some height above the 
disk. The result that the solutions follow the stable trends, 
confirms intuitive expectation. 

However, the existence of the critical magnetic diffusiv- 
ity, now illustrated clearly for the first time in comparison 
with the analytical solution of the closely related ideal-MHD 
problem, sets a limit for such intuitive reasoning. 



A general conclusion for the resistive MHD simulations 
of jets is that they are similar to the ideal MHD solutions, 
for a finite range of the parameter of magnetic diffusivity. In 
this range, they reach a well defined stationary state. This 
also extends to the numerical resistivity implicitly present in 
codes. In this respect our result confirms that in numerical 
simulations with reasonable resolution, the result should not 
differ significantly from the ideal MHD solution. Departure 
of the solution from the ideal-MHD regime seems to occur, 
at least for simple, smooth initial and boundary conditions, 
only for larger values of the magnetic diffusivity, a few orders 
of magnitude above the level of numerical magnetic diffusiv- 
ity. This regime of our solutions will be investigated in more 
detail in a following study. 

Here we presented an application of our results to the 
case of jets associated with young stellar objects. Evidently, 
these results are scalable to various other astrophysical 
cases, as well. 
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Figure 12. Values of i?^ for the final state with r\ = 0.03, at 
resolution R X Z = 256 X 512. The top panel shows contours of 
log 10 Rp (six contour lines per decade). It is also shown in colour 
grading. The bottom panel shows values of Rp along the vertical 
line indicated in the top panel. We see that the value just above 
the disk is Rp = 8.0 in the slice taken here. In the whole domain, 
the minimum value is Rp = 4. 
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